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Abstract 

A  new  approach  to  the  problem  of  computing  dense  optical  flow  fields  in  an  image 
sequence  is  presented.  Standard  formulations  of  this  problem  require  the  computation¬ 
ally  intensive  solution  of  an  eUiptic  partial  differential  equation  which  arises  from  the 
often  used  “smoothness  constraint”  type  regularization.  We  utilize  the  interpretation 
of  the  smoothness  constraint  as  a  “fractal  prior”  to  motivate  regularization  based  on 
a  recently  introduced  class  of  multiscale  stochastic  models.  The  solution  of  the  new 
problem  formulation  is  computed  with  an  efficient  multiscale  algorithm.  Experiments 
on  several  image  sequences  demonstrate  the  substantial  computational  savings  that  can 
be  achieved  due  to  the  fact  that  the  algorithm  is  non-iterative  and  in  fact  has  a  per 
pixel  computational  complexity  which  is  independent  of  image  size.  The  new  approach 
also  has  a  number  of  other  important  advantages.  Specifically,  multiresolution  flow  field 
estimates  are  available,  allowing  great  flexibility  in  dealing  with  the  tradeoff  between 
resolution  and  accuracy.  Error  covariance  information  is  also  available,  which  is  of  con¬ 
siderable  use  in  assessing  the  accuracy  of  the  estimates.  Finally,  the  algorithm  is  an 
excellent  “pre-conditioner”  for  the  iterative  algorithms  associated  with  the  smoothness 
constraint  problem  formulation.  Indeed,  its  usefulness  should  extend  to  a  wide  vari¬ 
ety  of  estimation  problems  including  many  involving  Gauss  Markov  random  fields  and 
spatial  processes  defined  as  the  solution  of  noise  driven  partial  differential  equations. 
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1  Introduction 


The  apparent  motion  of  brightness  patterns  in  an  image  is  referred  to  as  the  optical  flow 
[14].  In  computationed  vision,  opticed  flow  is  an  important  input  into  higher  level  vision  al¬ 
gorithms  performing  tasks  such  as  segmentation,  tracking,  object  detection,  robot  giiidance 
and  recovery  of  shape  information  [1,  19,  22,  25,  29].  In  addition,  methods  for  computing 
optical  flow  are  an  essential  peirt  of  motion  compensated  coding  schemes  [3,  33]. 

In  this  paper,  we  present  a  new  approach  to  the  problem  of  computing  optical  flow. 
Standard  formulations  of  this  problem  require  the  computationedly  intensive  solution  of  an 
elliptic  partial  diflerential  equation  which  arises  from  the  often  used  “smoothness  constraint” 
regularization  term.  We  utilize  the  interpretation  of  the  smoothness  constraint  as  a  “fractal 
prior”  to  motivate  regularization  based  on  a  recently  introduced  class  of  multiscale  stochastic 
models.  These  models  are  associated  with  efficient  multiscale  smoothing  algorithms,  and 
experiments  on  several  image  sequences  demonstrate  the  substantial  computational  savings 
that  can  be  achieved  through  their  use. 

Our  approach  is  most  easily  understood  by  comparing  it  to  the  optical  flow  problem 
formulation  originally  proposed  by  Horn  and  Schunck  [14].  As  they  discuss,  information 
about  the  optical  flow  field  can  be  obtained  from  the  image  sequence  by  making  the  as¬ 
sumption  that  changes  in  image  brightness  are  due  only  to  motion.  This  leads  to  the  so 
cedled  brightness  constraint  equation  [14]: 

d  S 

0  =  Z2,  t)  =  -^Eizi,  Z2,  t)  -I-  VE{zu  Z2-,  t)  •  ®(2l,  22,  t)  (1) 

where  E{zi,  22,  t)  is  the  image  intensity  as  a  function  of  time  t  and  space  (21, 22),  ®(2i,  22,  t) 
is  the  optical  flow  vector  field,  and: 


The  brightness  constraint  equation  (1)  does  not  completely  specify  the  flow  field  since 
it  provides  oidy  one  linear  constraint  for  the  two  unknowns  at  each  point.  This  is  usually 
referred  to  as  the  aperture  problem  [14].  One  way  to  obtain  a  unique  solution  is  to  regularize 
the  problem  by  imposing  an  additiontd  smoothness  constraint.  Specificeilly,  one  formulates 
the  following  optimization  problem  [14]: 

isc  =argmin  j  j  R~^  dz\dz2  (4) 

The  smoothness  constraint  is  captured  by  the  second  term  which  penalizes  large  gradients 
in  the  optical  flow.  The  constant  R  allows  one  to  tradeoff  between  the  relative  importance 
in  the  cost  function  of  the  brightness  and  smoothness  constraint  terms.  We  refer  to  the 
optical  flow  estimate  obtained  from  (4)  as  the  smoothness  constraint  (SC)  solution  to  the 
problem  of  computing  optical  flow. 

One  of  the  major  problems  associated  with  the  formulation  in  (4)  is  that  it  leads  to 
computationally  intensive  algorithms.  Specifically,  one  can  show  that  the  solution  of  (4) 
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Figrire  1:  Depiction  of  three  fields  which  are  equally  favored  by  the  smoothness  constraint, 
illustrating  how  this  penalty  provides  a  fractal  prior  model  for  the  optical  flow. 

satisfies  an  elliptic  partial  differential  equation  (PDE)  [14].  Discretization  of  this  PDE  leads 
to  a  sparse  but  extremely  large  set  of  linear  equations.  The  equations  are  typically  solved 
using  iterative  approaches,  which  require  increasing  numbers  of  iterations  as  the  image  size 
grows.  One  of  the  first  iterative  approaches  used  was  the  Gauss-Seidel  relaxation  algorithm 
[14,  27]  which  is  extremely  simple,  but  converges  very  slowly.  Terzopoulos  [31]  proposed 
the  use  of  multigrid  approaches  and  reported  a  factor  of  7  reduction  in  computation  over 
the  Gauss-Seidel  approach.  Successive  over-relaxation  (SOR)  algorithms  [16]  also  provide 
significant  computational  improvement  over  GS  approaches  and  have  been  successfully  used 
in  [21,  23,  24]. 

The  meiin  purpose  of  this  paper  is  to  address  the  computational  issue  discussed  above. 
To  do  this,  we  need  to  analyze  the  smoothness  constraint  in  more  detail.  Note  that  the 
penalty  associated  with  the  smoothness  constraint  term  in  (4)  is  equal  to  the  integral 
of  the  squared  norm  of  the  field  gradient  over  the  image  plane.  In  a  one- dimensioned 
context,  such  a  constraint  would  penedize  each  of  the  (one-dimensional)  fields  in  Figure  1 
equally.  Intuitively,  the  smoothness  constraint  has  a  fractal  nature,  and  in  fact  this  can 
be  demonstrated  in  a  much  more  precise  sense.  As  discussed  in  Rougee  et.  al.  [23,  24], 
the  optical  flow  problem  formulation  in  (4)  has  an  equivalent  formulation  in  an  estimation- 
theoretic  context.  Specifically,  let  us  make  the  assumption  that  the  two  components  of  the 
optical  flow  field  are  independent,  two-dimensional  Brownian  motions^.  Next,  suppose  that 

*More  precisely,  since  we  do  not  want  to  bias  the  optical  flow  estimates  towards  zero,  we  only  assume 
that  the  gradients  of  the  optical  flow  field  components  are  equal  to  the  gradients  of  the  Brownian  motion 
processes.  This  avoids  placing  a  constraint  on  the  DC  (i.e.  average)  value  of  the  optical  flow  and  focuses 
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(5) 


we  interpret  the  brightness  constraint  as  a  noisy  measurement  of  flow  field: 

-  -^2,  t)  =  VE{zi,  Z2,  t)  •  x{zi,  Z2,  t)  +  u(zi,  22,  0 

where  w(zi,  22,  t)  is  zero  mew  white  Gaussian  noise  with  intensity  R.  Then,  the  statistically 
optimal  estimate  of  the  flow  field,  given  the  measurements  (5)  and  the  Browtiian  motion 
prior  model,  is  the  same  as  the  optical  flow  estimate  given  by  (4).  Thus,  the  smoothness  con¬ 
straint  and  the  estimation- theoretic  formulations  are  eqxii valent.  The  estimation-theoretic 
interpretation  simply  allows  us  to  interpret  the  smoothness  constraint  as  a  Browniw  mo¬ 
tion  model.  In  one-dimension,  Brownian  motion  is  a  statistically  self-similar,  fractal  process 
with  &IIP  generalized  spectrum  [18],  wd  for  this  reason  the  smoothness  constraint  is  often 
referred  to  as  a  “fractal  prior”  [30]. 

Given  that  the  smoothness  constraint  prior  model  has  been  introduced  solely  for  regu¬ 
larization  purposes,  we  are  led  to  the  idea  of  replacing  it  with  another  prior  model  which  is 
simileir  in  its  fractal  nature,  but  which  leads  to  a  computationally  more  attractive  problem 
formtdation.  Replacing  the  prior  model  simply  means  that  we  want  to  change  the  regu¬ 
larization  term  in  (4).  The  interpretation  of  the  smoothness  constraint  as  a  fractal  prior 
suggests  that  other  prior  models  which  also  have  self-siznilar  characteristics  would  provide 
comparable  flow  estimates.  In  this  paper,  we  demonstrate  an  approach  which  substitutes 
a  fractal-like  class  of  prior  models  recently  introduced  in  [7,  8,  9,  11]  for  the  smoothness 
constraint  term.  The  structure  of  the  new  penalty  term  leads  to  an  extremely  efficient 
algorithm  for  the  computation  of  optical  flow  estimates.  The  algorithm  is  not  iterative 
and  in  fact  requires  a  fixed  number  of  floating  point  operations  per  pixel  independent  of 
image  size.  Thus,  since  all  of  the  known  methods  for  solving  the  smoothness  constrednt 
problem  formulation  have  per  pixel  computationad  complexities  that  grow  with  image  size, 
the  computational  savings  associated  with  the  new  approach  increases  as  the  image  size 
grows. 

In  addition  to  the  computational  savings,  there  are  several  other  advwtages  that  the 
new  approach  has.  First,  the  models  provide  a  representation  of  the  optical  flow  field  at 
multiple  resolutions.  For  this  reason,  we  say  that  they  provide  a  multiscale  regularization 
(MR)  of  the  optical  flow  problem,  and  we  refer  below  to  the  MR  eilgorithm  and  solution. 
The  multiscale  representation  of  the  flow  field  should  provide  great  flexibility  in  dealing 
with  the  tradeoff  between  accuracy  and  resolution.  Specifically,  one  can  expect  to  obtain 
higher  accuracy  at  lower  resolutions,  md  one  should  be  able  to  use  the  flow  estimates  at  one 
resolution  to  decide  if  computational  resources  should  be  used  to  compute  representations 
at  finer  resolutions. 

Second,  because  we  develop  the  approach  in  an  estimation-theoretic  context,  the  flow 
estimates  have  associated  error  covariance  matrices  which  provide  information  about  their 
quality.  This  information  would  be  essential  to  addressing  the  tradeoff  between  resolution 
and  accuracy  as  discussed  above,  and  may  also  be  useful  to  higher  level  vision  algorithms 
which  need  to  combine  information  in  a  rational  way  from  a  variety  of  sources  [26] 

A  third  advantage  of  the  approach  is  that  it  is  an  excellent  pre-conditioner  for  algorithms 
which  compute  the  smoothness  constraint  solution.  We  stress  that  the  optical  flow  estimates 
obtained  with  the  new  problem  formulation  will  not  exactly  equal  those  given  by  (4),  since 

only  on  imposing  a  preference  for  smoothness  in  the  flow. 
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the  prior  model  is  not  exactly  the  same.  Experimental  evidence  in  Section  3  shows  that  the 
difference  between  the  SC  and  MR  flow  estimates  consists  of  mostly  high  spatial  frequency 
components.  These  are  precisely  the  components  which  can  be  quickly  removed  by  the 
iterative  algorithms  computing  a  smoothness  constraint  solution.  Thus,  if  one  requires  the 
estimate  based  on  the  smoothness  constraint,  it  can  be  obtained  by  using,  for  instance,  a 
“pre-conditioned  SOR”  algorithm,  which  utilizes  the  MR  solution  as  an  initicd  estimate  of 
the  optical  flow,  amd  we  present  restdts  that  demonstrate  the  efliciency  of  this  approach. 

Finally,  our  approach  should  be  applicable  in  a  substantially  more  general  setting.  In 
particidar,  it  may  provide  a  computationally  attractive  alternative  to  standard  approaches 
to  the  broad  class  of  estimation  problems  in  which  the  underlying  process  to  be  estimated 
is  modeled  as  a  Gaussian  Markov  random  field  or  is  obtained  as  the  solution  of  noise  driven 
partial  differential  equations,  or  in  which  a  “smoothness  constraint”  type  regularization  is 
employed. 

This  paper  is  organized  as  follows.  In  Section  2  we  discuss  in  more  detail  an  estimation- 
theoretic  interpretation  of  the  optical  flow  formulation  in  (4)  and  develop  our  new  approach 
to  the  computation  of  optical  flow.  Section  3  presents  experimental  results  on  several  real 
and  synthetic  image  sequences.  Section  4  provides  further  discussion  sind  conclusions. 

2  Multiscale  Regularization 

In  the  first  part  of  this  section  we  develop  a  discrete  formulation  of  the  optical  flow  problem, 
and  discuss  in  more  detcdl  the  estimation- theoretic  interpretation  of  it.  We  then  illustrate 
precisely  how  the  smoothness  constraint  can  be  interpreted  as  a  prior  model  for  the  flow  field, 
jind  how  it  can  be  replaced  by  another,  similar  prior  model  which  leads  to  a  computationally 
more  attractive  problem  formulation.  The  general  class  of  prior  models  we  use  is  then 
introduced  along  with  an  algorithm  for  finding  the  solution  of  the  new  optical  flow  problem 
formtilation. 

2.1  An  Estimation-Theoretic  Interpretation  of  the  Optical  Flow  Problem 

We  start  by  introducing  the  following  notation.  Define: 

d 

y{zi,z2)  =  -—E{zi,Z2,t)  (6) 

C{zi,Z2)  =  VE{zuZ2,t)  (7) 

The  brightness  constraint  equation  (1)  can  then  be  written: 

yizi,Z2)  =  C{Z1,Z2)'X{Z1,Z2)  (8) 

where  the  time  dependence  of  the  equations  has  been  suppressed.  In  practice,  brightness 
measurements  are  only  available  over  a  discrete  set  of  points  in  space  and  time.  Thus,  the 
temporal  and  spatial  derivative  terms  in  the  brightness  constraint  equation  (8)  must  be 
approximated  with  finite  differences,  and  the  opticed  flow  is  only  estimated  on  a  discrete 
space-time  grid.  There  are  a  number  of  important  issues  which  arise  due  to  the  discretiza¬ 
tion;  we  refer  the  reader  to  [6]  for  a  detailed  discussion.  We  will  assume  here  that  the  optical 
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flow  is  to  be  estimated  on  the  set  22)1^1  =  ih,  Z2  =  jh\  i,j  £  {!,•••,  2^}}  where  h  is 
the  grid  spacing  and  M  is  an  integer.  The  asstimption  that  the  lattice  is  square  and  that 
the  number  of  rows  is  equal  to  a  power  of  two  makes  the  subsequent  development  easier,  but 
this  is  not  essential  as  we  discuss  in  Appendix  A.  Abusing  notation,  we  write  the  brightness 
constraint  (8)  at  the  point  {ih,jh)  as: 

yii,j)  =  C{i,j)-x{i,j)  (9) 

where  y{i,j)  is  the  measmed  temporal  brightness  derivative,  x(t,j)  is  the  optical  flow,  and 
C{i,j)  is  the  spatial  gradient  of  the  image  brightness  at  grid  point  {ih,jh). 

The  brightness  constraints  (9)  at  all  grid  points  can  be  grouped  into  one  large  set  of 
linear  equations  to  capture  the  optical  flow  information  contained  in  the  image  sequence. 
Deflning  x  as  the  vector  of  optical  flow  vectors  at  all  grid  points  (using,  say,  a  lexicographic 
ordering),  C  as  the  matrix  containing  the  corresponding  spatial  gradient  terms  C{i,j),  and 
y  as  the  vector  of  temporal  gradients,  we  can  write: 

y  =  Cx  (10) 


Then,  the  discrete  counterpart  of  (4)  is: 

±sc  =  argmin  ||y  -  Cx|||t_i  +  ||Lx||| 

X 

=  argmin  (y  —  Cx)^R“^(y  —  Cx)  +  x^L^Lx  (11) 

X 

where  the  matrix  L  is  a  discrete  approximation  of  the  gradient  operator  in  (4)  and  R  =  RI. 
The  regularization  term  x^L^Lx  makes  the  optimization  problem  (11)  well-posed.  In 
particular,  the  solution  of  (15)  satisfies  the  so  called  normal  equations  [28]: 

(C^R-^C-fL^L)xsc  =  (12) 

and  the  invertibility  of  (C^R~^C  -t-  L^L)  guarantees  that  xsc  is  unique.  The  normal 
equations  (12)  are  the  discrete  counterpart  of  the  partial  differential  equation  that  arises 
from  (4). 

An  estimation-theoretic  formulation  of  the  optimization  problem  in  (11)  can  now  be 
developed,  aind  we  will  use  it  to  show  that  the  statistically  optimal  estimate  of  the  optical 
flow,  given  a  particular  set  of  measurements,  is  identical  to  the  smoothness  constrciint 
solution  given  in  (11).  Specifically,  given  the  measurements: 

y  =  Cx-fv  (13) 

0  =  Lx  +  w  (14) 

and  the  assumptions  that®  v  ~  A/’(0,  R)  and  w  ~  M{0,I),  the  measurement  vector  y  = 
[y^|0]^  is  conditionally  Gaussian,  and  the  maximum  likelihood  estimate  [32]  of  x  is: 

XML  =  argmax  p(y|x) 

X 

*The  notation  z  ~  A/'(m,  A)  means  that  z  has  a  Gaussian  distribution,  with  mean  m  and  variance  A. 
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=  argmin  —  logp(y|x) 


argmin  [constant]  +  — 

X  2 


y 

0 


c 

L 


T  r 


R 


argmin  (y  -  Cx)^R  ^(y  -  Cx)  +  x^L^Lx 

X 

XSC 


(15) 


Thus,  the  msiximum  likelihood  problem  formulation  results  in  the  same  solution  as  the 
smoothness  constrsunt  formulation  when  L  is  used  to  define  an  additional  set  of  noisy 
measurements.  The  main  point  here  is  that  by  formulating  the  problem  in  this  estimation- 
theoretic  framework,  we  can  use  (14)  to  interpret  the  smoothness  constreunt  as  a  prior 
probabilistic  model  for  the  flow  field.  Specifically,  we  can  rewrite  (14)  as: 


Lx  =  -w  (16) 

Recalling  that  L  is  an  approximation  to  the  gradient  operator,  we  see  that  (16)  is  nothing 
more  than  a  spatial  difference  equation  model  for  x  driven  by  the  spatial  white  noise  field 
w. 

To  some  extent  the  precise  form  of  this  prior  model  is  arbitrary,  and  thus  we  me  led  to 
the  idea  of  introducing  a  new  prior  model  which  is  similar  in  nature,  but  which  leads  to 
a  computationally  more  attractive  problem  formulation.  That  is,  we  want  to  change  the 
smoothness  constraint  term  x^L^Lx  in  (15)  to  something  similar,  say,  x^Sx  «  x^L^Lx 
(where  S  is  a  symmetric  positive  semi-definite  matrix)  such  that  the  resulting  optimization 
problem  is  easy  to  solve.  If  we  factor  S  as  S  =  L^L  then  we  can  interpret  the  new  constraint 
as  a  prior  probabilistic  model  just  as  we  did  with  the  smoothness  constraint.  In  addition, 
there  is  a  precise  interpretation  of  what  we  have  done  as  a  Bayesian  estimation  problem. 
Specifically,  if  S  is  invertible  and  we  let  A~^  =  S,  then  the  use  of  this  new  constraint  in 
place  of  the  smoothness  constraint  is  equivalent  to  modeling  the  flow  field  probabilistically 
as  X  ~  A),  since  in  this  case  the  Bayes’  least  squares  estimate  of  the  flow  field  x,  given 

the  measurements  in  (13)  is  given  by: 

^BLSE  =  arpmtn  (y  -  Cx)^R~^(y  -  Cx) -fx^A'^x  (17) 

X 

which  corresponds  to  (15)  with  a  different  prior  model  term.  The  normal  equations  corre¬ 
sponding  to  (17)  are  given  by: 

(C^R-^C -f-A-^  )xb£s£;  =  C^R-V  (18) 


Comparison  of  the  problem  formulations  (11)  and  (17),  or  of  the  normal  equations  (12) 
and  (18),  makes  it  apparent  how  the  two  problem  formulations  are  related.  Note  that  an 
analogous  Bayesian  interpretation  can  apparently  be  given  to  the  smoothness  constraint 
formulation  (11),  (12),  with  the  corresponding  prior  model  for  optical  flow  given  by  x  ~ 
^*(0,  (L^L)”^).  Recall,  however,  that  L  is  an  approximation  to  the  spatial  gradient  operator 
and  thus  is  not  invertible  since  operating  on  constants  with  this  operator  yields  zero.  The 
probabilistic  interpretation  of  this  is  that  the  model  (16)  places  probabilistic  constraints  on 
the  spatial  differences  of  the  optical  flow,  but  not  on  its  DC  value.  Indeed,  it  is  not  difficult 
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to  check  that  if  we  model  optical  flow  instead  as  x  ~  A/’(0,  (L^L  +  where  e  is  any 

Mbitrarily  small  positive  number,  then  L^L  +  el  is  indeed  invertible  and  the  DC  value  of  x 
has  a  prior  covariance  Po  on  the  order  of  1  /e,  so  that  Po  oo  as  e  0.  Thus,  the  original 
smoothness  constraint  formulation  in  essence  assumes  an  infinite  prior  covariance  on  the 
DC  value  of  optical  flow.  The  alternate  model  developed  in  the  next  section  has  a  similar 
parameter,  Po,  representing  the  DC  variance,  which  can  similarly  be  set  to  oo. 

The  choice  of  the  new  prior  model  is  now  clearly  at  the  heart  of  the  problem.  Recalling 
that  the  smoothness  constraint  has  the  interpretation  as  a  “fractal  prior”,  we  would  like  to 
choose  a  prior  model  which  also  has  fractal-like  chsiracteristics.  A  natural  way  to  specify 
such  models  is  to  explicitly  represent  the  optical  flow  field  at  multiple  scales.  A  stochastic 
modeling  framework  which  allows  us  to  do  this,  and  which  also  leads  to  efficient  algorithms 
for  solving  (17),  (18),  is  described  in  the  next  section. 

2.2  A  Class  of  Multiscale  Models 

The  models  we  utUize  to  replace  the  smoothness  constraint  prior  model  were  recently  in¬ 
troduced  in  [7,  8,  9,  11].  The  models  represent  the  flow  field  at  multiple  scales,  i.e.  for  a 
set  of  scales  m  =  0, . . . ,  M,  with  m  =  0  being  the  coarsest  scale  and  m  =  M  the  finest 
scale,  we  define  a  set  of  optical  flow  fields  indexed  by  scale  and  space,  namely  Xm{i,j).  At 
the  scale,  the  field  consists  of  4”*  flow  vectors,  as  illustrated  in  Figure  2,  capturing 
features  of  the  optical  flow  field  discernible  at  that  scale  (i.e.  finer  resolution  features  of  the 
field  appear  only  in  finer  scale  representations).  Thus,  the  coarsest  version  of  the  flow  field 
consists  of  just  a  single  vector  corresponding  to  the  average  value  of  the  optical  flow  over  the 
entire  spatial  domain  of  interest,  and  successively  finer  versions  consist  of  a  geometrically 
increasing  number  of  vectors.  At  the  finest  level,  the  flow  field  is  represented  on  a  grid  with 
the  same  resolution  as  the  image  brightness  data.  In  particular,  XM{hj)  corresponds  to  the 
optical  flow  vector  x{i,j)  in  (9). 

Abstractly,  we  are  representing  the  flow  field  on  the  quadtree  structure  illustrated  in 
Figure  3.  Pyramidal  data  structures  such  as  the  quadtree  naturally  arise  in  image  process¬ 
ing  algorithms  which  have  a  multiresolution  component.  For  instance,  successive  filtering 
and  decimation  operations  lead  to  images  defined  on  such  a  hierarchy  of  grids  in  the  Lapla- 
cian  pyrjimid  coding  algorithm  of  Burt  £ind  Adelson  [5]  and  in  the  closely  related  wavelet 
transform  decomposition  of  images  [17].  Also,  the  multigrid  approaches  to  low  level  vision 
problems  discussed  by  Terzopoulos  [31]  involve  relaxation  on  a  similar  sequence  of  grids. 

The  model  we  introduce  in  this  section  describes  in  a  probabilistic  manner  how  the 
optical  flow  field  ®(i,  j)  —  xjaih  j)  is  constructed  by  adding  detail  from  one  scale  to  the  next. 
Just  as  the  smoothness  constraint  prior  model  (16)  describes  probabilistic  constraints  among 
values  of  the  optical  flow  at  different  spatial  locations,  our  multiscale  model  describes  such 
constraints  among  values  at  different  scales.  That  is,  our  model  describes  the  probabilistic 
evolution  of  Xm(i,j)  as  the  scale  m  evolves  from  coarse  to  fine.  For  notational  convenience 
in  describing  such  models,  we  denote  nodes  on  the  quadtree  with  a  single  abstract  index 
s  which  is  associated  with  the  3-tuple  {m,i,j)  where,  again,  m  is  the  scale  and  {i,j)  is 
a  spatial  location  in  the  grid  at  the  scale.  It  is  also  useful  to  define  an  upward  shift 
operator  7.  In  particular,  the  parent  of  node  s  is  denoted  37  (see  Figure  3).  We  note  that 
the  operator  7  is  not  one-to-one;  it  is  in  fact  four-to-one  since  each  node  wiU  have  four 
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Figure  2:  The  structure  of  a  multiscale  optical  flow  field  is  depicted.  The  components  of  the 
field  are  denoted  where  m  refers  to  the  scale  and  the  pair  {i,j)  denotes  a  particular 

grid  location  at  a  given  scale.  At  the  coeirsest  scale,  there  is  a  single  flow  vector  and,  more 
generally,  at  the  scale  there  are  4"*  vectors. 
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Root  node 


m  =  0 


m=  1 


m  =  2 


Figure  3:  Quadtree  structiure  on  which  the  multiscale  processes  are  defined.  The  abstract 
index  s  refers  to  a  node  in  the  quadtree;  sy  refers  to  the  parent  of  node  a. 


“offspring”  at  the  next  scale.  For  instance,  if  a  corresponds  to  any  of  the  nodes  in  the  upper 
left  quadrant  of  the  second  level  grid  (see  Figure  2),  i.e.  nodes  (2, 1, 1),  (2, 2, 1),  (2, 1,2)  or 
(2, 2, 2),  then  sy  corresponds  to  their  parent  on  the  first  level,  namely  node  (1, 1, 1). 

We  are  now  in  a  position  to  describe  the  class  of  multiscale  models  which  describe  the 
evolution  of  multiscale  stochastic  processes  indexed  by  nodes  on  the  quadtree.  Specifically, 
a  stochastic  quadtree  process  a;(a)  is  described  recursively  by: 


iB(s)  =  ^(a)«(s7)  -h  B{s)w{s) 

(19) 

imder  the  following  assumptions: 

*0  ~  -^(0,  Po) 

(20) 

U7(s)  ~  ^/’(O,/) 

(21) 

The  vectors  x  and  w  are  referred  to  as  the  state  and  driving  noise  terms.  The  state  vari¬ 
able  *0  at  the  root  node  of  the  tree  provides  an  initial  condition  for  the  recursion.  The 
driving  noise  is  white  in  both  space  and  scede,  and  is  uncorrelated  with  the  initial  con¬ 
dition.  Interpreting  each  level  as  a  representation  of  a  two-dimensional  field,  we  see  that 
(19)  describes  the  evolution  of  the  process  from  coarse  to  fine  scales.  The  term  A(s)x(sy) 
represents  interpolation  down  to  the  next  level,  and  B(s)w(s)  represents  higher  resolution 
detail  added  as  the  process  evolves  from  one  scale  to  the  next.  In  the  application  of  interest 
here,  a:(3)  =  Xtn(i,j),  where  s  =  (Tn,i,j),  and  thus  A,B  E  Such  a  model  corresponds 
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in  essence  to  a  first-order  recursion  in  scale  for  optical  fiow.'* 

Measurements  of  the  finest  level  optical  flow  field  are  available  from  the  brightness 
constraint.  In  particidar,  at  a  particidar  point  {i,j)  at  the  finest  level  M,  we  have  a 
measurement  equation  corresponding  to  that  in  (9): 

y{hj)  =  C{iJ)xM{hj)  +  v{i,j)  (22) 

v{i,j)  ~  A/'(0,i2)  (23) 

where  C7(i,j)  €  and  the  white  Gaussian  observation  noise  is  assumed  to  be  independent 
of  the  initial  condition  xo  and  the  driving  noise  w  in  (19)  -  (21).  Of  course,  we  can  group 
the  state  variables  a:(s)  at  the  finest  level  into  a  vector  xm  as  well  as  the  corresponding 
measurements  y(s)  and  spatial  gradient  terms  (^(s)  in  the  same  way  as  we  did  to  get  (10): 

y  =  Cxjif  ■+•  V  (24) 

V  ~  A/'(0,R)  (26) 

We  now  have  exactly  the  framework  which  led  to  the  statement  of  (17)  as  an  alternative 
to  the  smoothness  constraint  formulation  (15).  In  particular,  the  modeling  equations  (19) 
-  (21)  indicate  that  at  the  finest  level  of  the  quadtree,  the  flow  field  vectors  will  be  a  set 
of  jointly  Gaussian  random  variables  xm  ~  A/’(0,  A),  where  A  is  implicitly  given  by  the 
parameters  in  (19)  -  (21),  and  a  set  of  noisy  measurements  given  by  (24).  The  Bayes’  least 
squares  estimate  of  xm  given  the  measurements  in  (24)  and  the  prior  model  is: 

xjwr  =  argmin  (y  -  Cxjjf)^R~^(y  -  Cxm)  +  x|^A"^xm  (26) 

The  multiscale  modeling  framework  thus  provides  an  alternative  to  the  smoothness  con¬ 
straint  formulation  of  (11)  or  (15). 

What  remains  to  be  done  are  (1)  to  specify  a  model  within  this  class  that  has  charac¬ 
teristics  similar  to  those  of  the  smoothness  constraint  prior  model,  and  (2)  to  demonstrate 
why  the  use  of  this  alternate  multiresolution  formtdation  is  of  any  interest.  We  defer  the 
latter  of  these  to  the  next  section  and  focus  here  on  the  former.  In  particular,  for  our 
midtiscale  model  based  on  (19)  -  (21)  to  approximate  the  smoothness  constraint  prior  we 
would  like  to  choose  our  model  parameters  so  that  we  have  L^L  «  A“^.  The  observation 
that  the  prior  model  implied  by  the  operator  L  in  (15)  corresponds  to  a  Brownian  motion 
“fractal  prior”  suggests  one  approach  to  choosing  the  model  parameters.  In  particular, 
the  one-dimensional  Brownian  motion  has  a  1//^  generalized  spectrum  [18].  It  has  been 
demonstrated  that  such  processes  are  well  approximated  by  multiscale  models  such  as  ours 
in  one  dimension  if  geometrically  decreasing  powers  of  noise  are  added  at  each  level  m  of  the 
process  [9,  34].  In  particular,  this  motivates  the  choice  of  ^(s)  =  64“  *  I  in  (19),  where  I 

is  the  2x2  identity  matrix,  and  where  6  and  ft  are  scalar  constants.  The  constant  6  directly 
controls  the  overall  noise  power  in  the  process.  Also,  as  discussed  in  [34],  the  choice  of  ft 

*More  generally,  higher-order  recursions  in  scale  can  be  captured,  just  as  in  standard  state  space  models, 
by  increasing  the  order  of  the  model,  i.e.  the  dimension  of  ®(s).  In  this  case  the  actual  optical  How  at  node 
3  would  correspond  to  a  subset  of  the  components  of  a:(a),  with  the  remainder  of  *(«)  devoted  to  capturing 
the  memory  in  the  multiscale  recursion.  In  this  paper,  however,  we  restrict  ourselves  to  the  simple  first  order 
recursion. 
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controls  the  power  law  dependence  of  the  generalized  spectrum  of  the  process  at  the  finest 
resolution  as  well  as  the  fractal  dimension  of  its  sample  paths.  Specifically,  this  spectrum 
has  a  1//''  dependence.  Thus,  the  choice  of  ju  =  2  would  correspond  to  a  Brownian-like 
fractal  process.  To  achieve  greater  flexibility  in  both  the  modeling  and  estimation,  we  allow 
/X  to  be  a  peireoneter  that  can  be  varied.  In  addition,  receiU  that  in  the  smoothness  con¬ 
straint  formulation,  L^L  was  not  invertible  because  of  the  implicit  assumption  of  infinite 
prior  variance  on  the  DC  value  of  the  optical  flow  field.  In  our  multiscale  regularization 
context,  this  would  correspond  to  setting  Pq  equal  to  infinity  in  (20).  This  can  be  done 
without  diffictilty  in  the  estimation  algorithms  described  next,  but  we  have  found  that  it  is 
generally  sufficient  to  simply  choose  Pq  to  be  a  large  multiple  of  the  identity. 

In  closing  this  section  let  us  illustrate  some  of  the  types  of  sample  paths  that  result  from 
scalar  versions  of  multiresolution  models  of  this  type  (i.e.  where  a;(s)  is  a  sc2dar).  Sample 
paths  of  the  finest  level  of  a:(s)  from  several  scalar  processes  for  different  parameter  choices 
of  such  a  model  are  illustrated  in  Figures  4  to  6.  The  processes  are  shown  as  mesh  plots 
to  better  illustrate  the  changes  induced  by  choosing  different  model  parjimeters.  As  noted 
above,  the  parameter  A(s)  reflects  interpolation  from  coarse  to  fine  scales.  As  A(s)  goes 
to  zero,  the  scale-to-scale  “memory”  of  the  process  decreases,  and  in  fact  if  A(3)  =  0,  the 
process  will  be  white  noise  at  each  level.  The  effect  of  changing  A(s)  is  illustrated  in  Figures 
4  and  5.  The  parameter  in  the  driving  noise  term  also  effects  the  correlation  structure 
of  the  process.  As  ft  is  increased,  the  driving  noise  variance  is  reduced  from  scale-to-scede 
at  a  faster  rate.  Thus,  there  is  greater  interscale  correlation,  and  we  expect  the  process  to 
appear  more  structured.  The  effect  of  changing  ft  is  illustrated  in  comparing  Figures  4  and 
6.  As  these  figures  indicate,  the  choice  of  A(s)  =  1  leads  to  the  fine  scale  process  having 
obvious  memory  of  coeirse  scale  features,  while  the  larger  veJue  of  ft  in  Figure  6  leads  to 
sample  paths  with  smaller  fine  scJile  fluctuations  than  Figure  4. 

2.3  The  Multiscale  Regularization  Algorithm 

We  have  now  specified  a  class  of  models  which  will  eJlow  us  to  approximate  the  smooth¬ 
ness  constraint  prior  model.  The  simple  multiscale  structure  of  these  models  leads  to  very 
efficient  algorithms  for  computing  the  optimal  estimate  of  the  state  given  a  set  of  measure¬ 
ments.  One  of  these  edgorithms,  which  we  refer  to  as  the  Mriltiscale  Regularization  (MR) 
edgorithm,  was  developed  in  [7,  8,  9,  10]  for  one- dimensional  signals,  and  its  extension  to 
images  is  described  here. 

The  MR  algorithm  computes  the  Bayes  least  squMes  estimate  of  the  state  vectors  (19) 
given  the  measurements  (22)  in  two  steps.  The  first  step  is  an  upward  or  fine-to-coarse 
sweep  on  the  quadtree,  which  propagates  the  measurement  information  in  parallel,  level 
by  level,  from  the  fine  scale  nodes  up  to  the  root  node.  The  second  step  is  a  downward 
or  coarse-to-fine  sweep  which  propagates  the  measurement  information  back  down,  and 
throughout  the  tree.  The  result  is  the  least  squares  estimate  of  the  state  x(a)  at  each  node 
based  on  all  of  the  data.  The  details  of  the  upward  emd  downward  sweeps  are  given  below 
and  are  discussed  in  much  greater  detedl  in  [9,  10]. 

To  begin,  note  first  that  the  measurement  model  (22),  allowing  for  the  possibility  of 
spatially  varying  noise  intensity,  can  be  written  in  the  form: 

y(s)  -  C(s)x(s)  +  v(s)  (27) 
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Figure  6:  Finest  scale  of  a  scalar  quadtree  process  seimple  path:  -4(s)  =  1,5(3) 


»(a)  ~  AC(0,JJ(»))  (28) 

In  the  context  of  the  optic£d  flow  estimation  problem,  measurements  are  taken  only  on  the 
finest  level,  corresponding  to  C(s)  =  0  unless  s  is  a  node  at  the  finest  level  (i.e.  uidess 
3  =  However,  in  the  more  general  modeling  framework  discussed  in  [9,  10],  the 

measurements  may  be  available  at  any  node,  and  the  noise  variance  may  vary  with  node  as 
in  (28).  We  present  here  this  more  general  algorithm  in  which,  in  addition,  ic,  y  and  w  may 
be  of  arbitrary  dimension. 

The  model  given  by  (19)  -  (21),  (27)  -  (28)  is  a  downward  model  in  the  sense  that  the 
recursion  starts  from  the  root  node  and  propagates  down  the  quadtree  from  coarse-to-fine 
scales.  In  order  to  describe  the  upward  sweep  of  the  MR  algorithm,  we  need  a  corresponding 
upward  model.  This  upward  model  is  equivalent  to  the  downward  model  in  the  sense  that 
the  joint  second  order  statistics  of  the  states  1(3)  and  measurements  y{3)  are  the  same.  The 
upward  model  is  given  by®  [8,  9]: 

*(37)  =  F(3)«(3)  —  44~^(3)5(3)u;(s)  (29) 

y{s)  =  C'(s)a:(3)  +  u(3)  (30) 

where: 

F(3)  =  P,^A^{s)Pr^  (31) 

®We  use  E[x]  to  denote  the  expected  value  of  the  random  variable  ss  and  E[®|y]  to  denote  the  expected 
value  of  *  given  y. 
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i6(s)  =  u;(s)  -  E[u7(s)|®(a)] 

(32) 

E[t2>(s)t5^(5)]  =  I-B^P-^B{s) 

(33) 

=  Q(^) 

(34) 

and  where  P,  =  E[a!(j)®r(5)]  is  the  variance  of  the  state  at  node  s  and  evolves  according 
to  the  Lyapanov  equation: 


P,  =  ^(a)P,^^^(3)  + 

To  proceed  further  we  need  to  define  some  new  notation. 

=  s  or  s'  is  a  descendant  of  s}  (36) 

n-"  =  r.\M  (37) 

*(s'|s)  =  E[*(s')|r.]  (38) 

ic(s'|s+)  =  E[*(s')|y+]  (39) 

P(s'|s)  =  E[(*(s')  -  i(s'|s))(®(sO  -  i(s'|s))r]  (40) 

P(s'|s+)  =  E[(®(s')  -  ®(a'|s+))(®(s')  -  *(s'|s+))r]  (41) 

where  the  notation  F,  \  {s}  means  the  node  s  is  not  included  in  the  set  Y+.  The  upward 
sweep  of  the  MR  algorithm  begins  with  the  initialization  of  ®(s|s+)  and  the  corresponding 
error  covariance  P(s|s+)  at  the  finest  level.  Then,  the  estimate  update  is  computed  via: 

®(s|s)  =  ®(j»|«4-)  +  JS:(s)[p(s)  -  C'(s)®(s(s4-)]  (42) 

P(s|s)  =  [/-fi:(s)C'(s)]P(s|s+)  (43) 

K{s)  =  P(s|3+)C'^(s)F“^(s)  (44) 

F(3)  =  C'(s)P(s|a+)C'^(s)  +  P(s)  (45) 

Denote  the  offspring  of  ®(s)  as  x(sai),t  =  For  the  quadtree  model,  g  =  4.  The 

updated  estimates  are  predicted  back  up  to  the  next  level: 

®(s|sai)  =  P(sai)®(sai|sai)  (46) 

P(s|sai)  =  P(sai)P(sai|sai)F^(sai)  +  Q(saf)  (47) 

Q(saf)  =  A“^(sai)P(5ai)Q(sai)B^(sai)j4~^(sai)  (48) 

The  predicted  estimates  are  then  merged: 

®(3|5+)  =  P(s|s+)^P“^(s(5ai)®(5|saf)  (49) 

»=i 

P(s|s+)  =  [(l-g)p-i+^p-i(a|sai)]-'  (50) 

*=i 


The  upward  sweep  given  by  the  update,  predict  and  merge  equations  proceeds  recursively 
up  the  quadtree.  At  the  top  of  the  tree,  one  obtains  the  smoothed  estimate  of  the  root 
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node,  that  is,  an  estimate  based  on  ail  of  the  data.  The  estimate  and  its  error  covariance 


eire  given  by; 

®*(0)  =  ®(010)  (51) 

P'(0)  =  P(0|0)  (52) 

where  the  superscript  s  denotes  the  fact  that  these  are  smoothed  estimates,  that  is,  they  are 
based  on  all  of  the  available  data.  The  smoothed  estimate  and  associated  error  covariance 
at  the  root  node  provide  initialization  for  the  downward  sweep,  which  is  given  by: 

i*(s)  =  ®(s|a)  +  J(s)[*'(s^)  -  *(s7|s)]  (53) 

P*{s)  =  P(s|s)  +  J(^)[P*(57)-P(^7k)]^"’(^)  (54) 

J{s)  =  Pis\s)F^{s)P--^{s^\s)  (55) 


The  estimates  x*(s)  at  the  finest  level  of  the  quadtree  provide  the  solution  to  (26).  The  form 
of  the  algorithm  we  have  specified  here,  which  generalizes  standard  Kalman  filtering  and 
smoothing  algorithms  to  the  multiscale  context,  obviously  assumes  that  the  state  covariance 
P,  is  well  defined  and  finite,  and  it  is  not  difficult  to  see  from  (35)  that  this  will  be  the  case 
if  Po  is  finite.  There  is,  however,  an  alternate  form  of  this  algorithm  presented  in  [9,  10] 
which  generalizes  so-called  information  form  algorithms  for  standard  state  space  models 
and  which  propagates  inverses  of  covariances.  In  this  alternate  form  it  is  straightforward 
to  accommodate  the  setting  of  Pq  to  infinity  (which  corresponds  to  P^^  =  0),  and  we  refer 
the  reader  to  [9,  10]  for  details.  As  mentioned  previously,  we  have  foimd  that  setting  Po 
to  a  large  but  finite  multiple  of  the  identity,  and  then  using  (42)  -  (50),  (53)  -  (55),  yields 
excellent  results. 

3  Experimental  Results 

3.1  Specification  of  the  Multiscale  Model 

To  specify  the  MR  algorithm  completely  we  need  to  choose  the  parameters  in  (19)  -  (21), 
(27)  -  (28).  We  utilize  the  following  parameterization  of  the  model: 


®(«) 

=  *(57)  -f-  (64  s  )«j('S) 

(56) 

=  C'(s)®(s)  -1-  t;(s) 

(57) 

to(3) 

~  A/'(0,J) 

(58) 

vis) 

~  A'(0,P(s)) 

(59) 

®0 

~  A/’(0,  pi) 

(60) 

where  J  is  a  2  x  2  identity  matrix.  From  (56)  and  (58)  we  see  that  the  two  components  of 
the  optical  flow  field  are  modeled  as  independent  sets  of  random  variables,  and  that  each 
will  have  a  fractal-like  characteristic  due  to  the  choice  of  the  driving  noise  gain  B{s)  (as 
discussed  in  the  previous  section).  We  view  p  and  b  as  free  model  parsimeters  which  can 
be  varied  to  control  the  degree  and  type  of  regularization  in  much  the  same  way  that  the 
parameter  R  in  the  smoothness  constraint  formulation  (4)  is  used  to  tradeoff  between  the 
data  dependent  and  regularization  terms  in  the  optimization  functional. 
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As  discussed  previously,  the  measurements  ^(s)  and  measurement  matrix  (^(s)  come 
directly  from  the  image  temporeil  and  spatial  gradients,  which  are  available  at  the  finest  level 
of  the  quadtree.  In  the  experiments  described  below,  we  use  a  simple  two-image  difference 
to  approximate  the  temporal  gradient.  The  spatied  gradient  is  computed  by  smoothing  the 
image  with  a  3  X  3  Gaussian  kernel  followed  by  a  central  difference  approximation.  The 
additive  noise  variance  is  given  by  R{s).  We  have  found  empirically  that  the  choice  R{s)  = 
maa!(||C'(s)|p,  10)  worked  well.  This  choice  effectively  penalizes  large  spatial  gradients, 
which  are  likely  points  of  occlusion  where  the  brightness  constraint  equation  will  not  hold 
[26].  The  parameter  p  in  the  prior  covariance  of  the  root  node  was  set  to  p  =  100.  The 
distribution  (60)  on  the  root  node  effectively  says  that  we  are  modeling  the  optical  flow  field 
components  as  zero  mean  r^mdom  processes.  The  prior  covariance  reflects  our  confidence 
in  this  assumption.  Since  we  do  not  believe  that  any  prior  assumption  on  the  mean  of 
optical  flow  field  components  can  be  justified,  we  set  the  parameter  p  such  that  the  implied 
st€tndard  deviation  is  much  larger  than  the  sizes  of  the  flow  fields  we  expect  to  see. 

We  compwe  our  approach  computationally  cmd  visually  to  the  the  Gauss-Seidel  (GS) 
and  successive  over-relaxation  (SOR)  ^dgorithms,  which  can  be  used  to  compute  the  solu¬ 
tion  of  the  smoothness  constraint  formulation  given  by  (11)  or  (15).  The  details  of  these 
algorithms  can  be  found  in  Appendix  B.  Straightforward  analysis  shows  that  the  GS  and 
SOR  algorithms  require  14  and  18  floating  point  operations  (flops)  per  pixel  per  iteration 
respectively.  The  number  of  iterations  required  for  convergence  of  the  iterative  algorithms 
grows  with  image  size  [16].  For  reasonable  size  images  (say,  512  x  512),  SOR  may  require 
on  the  order  of  hundreds  of  iterations  to  converge,  so  that  the  total  computation  per  pixel 
can  be  on  the  order  of  10®  —  10'*  flops.  On  the  other  hand,  the  MR  algorithm  requires  76 
flops  per  pixel  (see  Appendix  C).  Note  that  the  MR  algorithm  is  not  iterative.  Thus,  the 
computational  gain  associated  with  the  MR  algorithm  can  be  on  the  order  of  one  to  two 
orders  of  magnitude. 

3.2  Rotation  Sequence 

The  first  example  is  a  sequence  of  Gaussian  images  modulated  by  a  spatial  sinewave.  Specif¬ 
ically,  the  first  frame  intensity  pattern  is  given  by: 

E{zi,Z2,ti)  =  sin(atan(zi  —  23, Z2  —  28))  exp(— iz'Z'^z) 
zi  -  23 

"  [  Z2  -  28 

„  ^  1000  0 
[  0  500 

where  atan(zi,Z2)  is  a  27r  arctangent  (atan(0,l)  =  0,  atan(l,0)  =  — tt),  h  =  1  and  M  =  6 
(i.e.  the  image  lattice  is  64  X  64,  cf.  the  discussion  about  discretization  at  the  beginning  of 
Section  2.1).  The  second  frame  is  equal  to  the  first,  rotated  by  1  degree  about  pixel  (23,28). 
The  first  frame  and  actual  optical  flow  are  illustrated  in  Figures  7  and  8. 

Figure  9  illustrates  the  flow  computed  using  the  MR  algorithm.  The  computed  flow 
reflects  the  rotational  nature  of  actual  flow  field,  with  error  at  the  boimdaries  induced  by 
the  quadtree  structure  of  the  prior  model.  The  degree  of  “blockiness”  is  determined  by  the 


(61) 

(62) 

(63) 
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parameters  b  and  fx.  Increasing  the  parameter  b  raises  the  level  of  uncertainty  in  the  prior 
model,  which  implies  that  the  MR  algorithm  will  tend  to  provide  less  regularization  (and 
hence  less  blockiness).  A  larger  value  of  fi  increases  the  prior  model  correlation  between 
any  two  given  states,  thereby  inducing  more  blockiness  in  the  flow  estimates. 

Figure  10  illustrates  the  smoothness  constraint  (SC)  flow  estimates  computed  using  the 
SOR  algorithm  (the  SOR  algorithm  was  initialized  with  identically  zero  flow  estimates). 
The  estimates  required  50  SOR  iterations  to  obtain,  representing  a  factor  of  50/4.2  =  11.9 
more  computation  than  the  MR  algorithm.  Figure  11  illustrates  the  root  mean  square^ 
(rms)  error  in  the  flow  estimates  as  a  function  of  iteration  for  the  SOR  and  GS  algorithms. 
As  expected,  the  SOR  algorithm  is  signiflcantly  faster  than  the  GS  algorithm  (they  will 
converge  to  the  same  result  since  they  are  solving  the  same  partial  differential  equation). 
The  rms  error  in  the  MR  flow  estimates  is  depicted  as  a  straight  line,  since  the  algorithm 
is  not  iterative. 

The  MR  and  SC  flow  estimates  are  not  identical  due  to  differences  in  the  prior  models. 
If  there  is  particular  interest  in  obtaining  the  SC  solution,  the  question  arises  of  using  the 
MR  solution  as  cm  initial  guess  for  the  iterative  algorithms  which  compute  the  SC  solution. 
Note  that  the  difference  between  the  SC  and  MR  flow  estimates  is  associated  with  the  non¬ 
smooth,  high  frequency  aspects  of  the  MR  flow  at  block  edges.  It  is  precisely  these  high 
frequency  components  that  are  quickly  removed  by  SOR  or  GS  algorithms  computing  the 
the  smoothness  constraint  solution  and  suggests  that  the  MR  algorithm  would  provide  an 
excellent  pre-conditioner  for  the  iterative  algorithms.  Figure  12  illustrates  the  optical  flow 
estimates  computed  using  a  pre-conditioned  SOR  algorithm.  The  estimates  correspond  to 
5  iterations  of  the  SOR  algorithm  initialized  with  the  MR  flow  estimates  in  Figure  9.  In 
this  exzunple,  the  total  computation  required  by  the  pre-conditioned  SOR  approach  is  a 
factor  of  50/(4.2  +  5)  =  5.4  less  than  that  required  by  the  stcmdard  SOR  algorithm  (i.e.  the 
algorithm  initialized  with  identically  zero  flow  estimates). 

Figure  13  illustrates  the  rms  difference  between  the  smoothness  constraint  solution  and 
the  intermediate  values  of  the  GS,  SOR  and  pre-conditioned  SOR  estimates  as  a  function 
of  iteration’^.  The  error  plot  for  the  pre-conditioned  SOR  algorithm  begins  at  4.2  iterations 
to  take  into  account  the  initial  computation  associated  with  the  MR  algorithm  (which 
equals  4.2  SOR  iterations).  The  figure  demonstrates  that  the  pre-conditioned  SOR  approach 
provides  a  substantially  less  computationally  intensive  approach  to  finding  the  SC  flow 
estimates  even  for  this  small  size  problem. 

3.3  Yosemite  Sequence 

The  second  example  is  a  synthetic  256  x  256  image  sequence  which  simulates  the  view  from 
a  small  plane  flying  through  the  Yosemite  Valley®.  The  first  image  in  the  sequence  is  shown 
in  Figure  14  along  with  the  actual  flow  field  in  Figure  15.  The  flow  computed  via  the 
MR  algorithm  is  shown  in  Figure  16  and  the  smoothness  constraint  solution  is  shown  in 
Figure  17.  The  smoothness  constraint  flow  estimates  required  250  SOR  iterations  in  this 

®The  rms  error  is  only  one  measure  of  flow  field  quality  and  will  not  be  appropriate  for  all  applications 
(e.g.  segmentation  or  coding). 

’^Tlie  smoothness  constraint  solution  is  approximated  as  the  SOR  algorithm  optical  flow  estimates  after 
500  iterations. 

*This  sequence  was  synthesized  by  Lyn  Quam  of  SRI  International. 


19 


Figure  7:  First  frame  of  the  “rotation”  sequence. 
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Figure  11:  Rms  Error  Comparison  of  MR,  SOR  and  Gauss- Seidel  (GS)  algorithm  flow 
estimates  for  the  rotation  sequence. 


//////' 
/  /  y  /  /  ^ 

y  /  y  /  /  /■ 

y  y  y  /  /  ^ 

/  /  y  /  /  f 

/  f  f  t  I  , 

t  t  t  !  \  \ 

f  f  f  t  \  s 

1  1  \  \  N  V 

\  \  \  \  V 

\  \  V  V. 


- 

\ 

\ 

\ 

\ 

\ 

N. 

\ 

\ 

\ 

\ 

\ 

" 

\ 

\ 

\ 

\ 

1 

1 

S 

\ 

\ 

\ 

1 

1 

1 

N 

\ 

\ 

1 

1 

J 

1 

1 

*  ' 

\ 

1 

1 

1 

/ 

1 

/ 

/ 

• 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

' 

✓ 

/ 

/ 

/ 

/ 

/ 

/ 

—  — 

✓ 

✓ 

y 

y 

✓ 

y 

/ 

/ 

—  ^ 

✓ 

y 

y 

y 

/ 

/ 

y 

y 

y 

y 

/ 

/ 

y 

y 

y 

y 

/ 

y 

y 

y 

/ 

/ 

/ 

X 

y 

/ 

/ 

/ 

/ 

/ 

/ 

y 

/ 

/ 

/ 

/ 

/ 

/ 

y 

/ 

/ 

/ 

/ 

/ 

/ 

Figure  12:  Pre-conditioned  SOR  algorithm  flow  estimates:  R  =  10^,  5  iterations.  The 
pre-conditioned  SOR  algorithm  is  initialized  with  the  MR  flow  estimates. 
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Figure  13:  Rms  Difference  Comparison  illustrates  how  the  pre-conditioned  SOR,  SOR  and 
GS  algorithms  converge  to  the  smoothness  constraint  solution.  The  plots  show  the  rms 
difference  between  the  smoothness  constraint  solution  and  the  estimates  as  a  function  of 
iteration.  All  will  eventually  converge,  but  the  pre-conditioned  SOR  algorithm  converges 
much  faster  than  SOR  or  GS. 

example,  representing  a  factor  of  60  more  computation  than  the  MR  estimates.  Note  the 
substantial  increase  over  the  previous  example  in  the  number  of  iterations  required  for  the 
SOR  algorithm  to  converge.  The  niimber  of  iterations  required  for  convergence  depends  on 
several  things,  including  the  parameter  R,  the  image  gradient  cheiracteristics  and  the  image 
size.  Theoretical  analysis  in  [16]  shows  that  the  SOR  algorithm  requires  on  the  order  of  N 
iterations  for  an  JV  x  iV  image.  Thus,  we  expect  substantially  more  computational  savings 
as  the  image  size  increases. 

An  rms  error  comparison  of  the  algorithms  is  shown  in  Figure  18.  Neither  of  the  es¬ 
timates  coincides  with  the  actual  optical  flow,  but  they  do  have  comparable  rms  error  as 
in  the  previous  example.  In  addition,  the  figure  illustrates  the  computational  advantage  of 
the  MR  algorithm.  In  particular,  the  SOR  algorithm  is  stiU  reducing  the  rms  error  in  its 
flow  estimates  after  300  iterations,  at  which  point  the  MR  algorithm  requires  a  factor  of 
300/4.2  =  71.4  less  computation. 

Again,  there  may  be  interest  in  obtaining  the  solution  to  the  smoothness  constraint 
problem  formulation.  Figure  19  depicts  the  pre-conditioned  SOR  results  (based  on  initial¬ 
ization  with  the  MR  flow  estimates  in  Figm-e  16)  after  20  iterations.  Visually,  there  is 
almost  no  difference  between  the  pre-conditioned  SOR  estimates  and  the  estimates  shown 
in  Figure  17.  The  computational  gain  associated  with  the  pre-conditioned  algorithm  is  10.3. 
Figure  20  illustrates  how  the  GS,  SOR  and  pre-conditioned  SOR  algorithms  converge  to 
the  smoothness  constraint  solution.  For  any  given  number  of  iterations,  the  pre-conditioned 
SOR  estimates  are  substantially  closer  to  the  final  solution  than  the  GS  or  SOR  estimates. 

This  image  sequence  contains  a  problem  often  encountered  in  real  images:  regions  of 
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Figure  14:  First  frame  of  “Yosemite”  sequence. 

constant  intensity.  The  problem  is  that  the  lack  of  gradient  information  in  that  region 
implies  that  the  optical  flow  is  not  well  defined.  The  smoothness  constraint  and  multiscale 
prior  models  provide  a  means  of  interpolating  out  into  these  regions.  The  result  of  this  is  ap¬ 
parent  in  the  top  portion  of  Figures  16  and  17.  An  advantage  of  the  MR  formulation  is  that 
it  accomplishes  this  extrapolation  at  an  appropriately  coarser,  eind  hence  computationally 
simpler,  scale. 

3.4  Moving  Vehicle  Sequence 

The  third  example  is  a  real®  512  x  512  image  sequence  and  depicts  the  view  from  a  car 
driving  down  a  road.  The  first  image  in  the  sequence  is  illustrated  in  Figure  21.  The  MR 
and  SOR  flow  estimates  are  shown  in  Figures  22  -  23.  In  this  example,  the  MR  solution 
required  a  factor  of  800/4.2  =  190  less  computation. 

Since  the  true  optical  flow  is  not  available  (as  it  was  in  the  previous  simulated  examples), 
an  alternate  performance  metric  is  needed.  In  particular,  we  wiU  use  a  reconstruction  error 
metric,  which  is  often  used  in  contexts  in  which  one  is  interested  in  using  opticad  flow  for 
motion  compensated  coding.  This  metric  measures  the  mean  square  difference  between 
the  current  image  in  a  sequence  and  an  estimate  of  it  based  on  the  computed  optical  flow, 
the  previous  image,  and  a  bilinear  interpolation  scheme  [20].  The  optical  flow  used  is  that 
associated  with  the  current  image.  Essentially,  one  estimates  the  brightness  at  any  given 
point  by  using  the  optical  flow  to  project  that  point  back  to  the  previous  image.  In  general, 
that  point  will  not  be  on  the  image  plane,  and  the  bilinear  interpolation  is  reqiured. 

®The  sequence  was  provided  by  Saab-Scania. 
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Figure  17:  SOR  algorithm  flow  estimates:  R  =  600*,  250  iterations.  The  SOR  algorithm 
required  a  factor  of  60  more  computation  in  this  example. 
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Figure  18:  Rms  Error  Comparison  of  MR,  SOR  and  Gauss-Seidel  (GS)  algorithm  flow 
estimates  for  the  Yosemite  sequence. 
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Figure  19:  Pre-conditioned  SOR  estimates:  R  =  500®,  20  iterations. 


Figure  20:  Rms  Difference  Comparison  illustrates  how  the  pre-conditioned  SOR,  SOR  and 
GS  algorithms  converge  to  the  smoothness  constraint  solution.  The  plots  show  the  rms 
difference  between  the  smoothness  constraint  solution  and  the  estimates  as  a  function  of 
iteration.  The  pre-conditioned  SOR  algorithm  is  substantially  closer  to  the  final  smoothness 
constraint  result  after  emy  given  number  of  iterations. 
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Figure  21:  First  frame  of  “Moving  Vehicle”  sequence. 


Figure  24  depicts  the  rms  reconstruction  error  based  on  the  MR,  SOR  and  GS  ap¬ 
proaches  as  a  function  of  iteration  and  illustrates  that  two  approaches  provide  comparable 
reconstruction  error  results,  just  as  they  provided  comparable  rms  error  results  in  the  previ¬ 
ous  examples.  In  addition,  it  illustrates  the  computational  advantage  of  the  MR  algorithm. 
Indeed,  the  SOR  algorithm  is  stiU  improving  its  optical  flow  estimates  (with  respect  to  this 
metric)  after  200  iterations,  at  which  point  the  MR  algorithm  requires  a  factor  of  200/4.2 
=  47.6  less  computation. 

Again,  we  illustrate  how  the  MR  algorithm  can  be  used  as  a  pre-conditioner.  Figure  25 
depicts  pre-conditioned  SOR  optical  flow  estimates  after  30  iterations.  The  pre-conditioned 
SOR  results  are  nearly  indistinguishable  from  the  standard  SOR  estimates  in  Figure  23  and 
represent  a  factor  of  23  reduction  in  computational  cost.  Figure  26  provides  a  view  of  how 
the  GS,  SOR  and  pre-conditioned  SOR  algorithms  converge  to  the  smoothness  constraint 
estimates  and  illustrates  that  the  pre-conditioned  SOR  algorithm  converges  significantly 
faster  than  the  other  algorithms. 

Note  finally  that  the  niimber  of  iterations  required  for  convergence  in  this  example  is 
greater  than  that  in  the  previous  two,  again  supporting  the  theoretical  results  in  [16].  Thus, 
one  can  expect  to  achieve  greater  computational  savings  through  the  MR  algorithm  as  the 
image  size  grows. 

4  Conclusions 

We  have  presented  a  new  approach  to  the  computation  of  optical  flow.  A  new  problem 
formulation  is  developed  which  utilizes  the  “fractal  prior”  interpretation  of  the  smoothness 
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Figure  22:  MR  algorithm  flow  estimates:  b  =  0.025,  ft  =  0.5. 
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Figure  23:  SOR  algorithm  flow  estimates:  R  =  500^,  800  iterations.  The  SOR  algorithm 
required  a  factor  of  190  more  computation  in  this  example. 
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Figure  24:  Reconstruction  Error  Comparison  of  MR,  SOR  and  Gauss-Seidel  (GS)  algorithm 
flow  estimates  for  the  Saab  sequence. 
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Figure  25:  Result  of  post-processing  MR  estimates  with  30  iterations  of  the  SOR  algorithm. 
Total  computational  effort  is  equivalent  to  34.2  iterations  of  the  SOR  algorithm. 
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an  initial  condition  for  the  iterative  algorithms  computing  a  smoothness  constraint  solution. 


constraint  to  motivate  regularization  based  on  a  recently  introduced  class  of  multiscale 
stochastic  models.  The  approach  has  the  advsmtage  that  the  solution  of  the  new  problem 
formulation  can  be  computed  very  fast,  in  comparison  to  the  solution  of  formulations  based 
on  a  smoothness  constraint.  In  fact,  the  computational  savings  realized  are  on  the  order 
of  a  factor  of  10  to  100.  Larger  gedns  are  associated  with  larger  images  since  the  iterative 
approaches  associated  with  the  smoothness  constraint  solution  take  longer  to  converge  as 
the  image  grows,  whereas  the  per  pixel  computation  associated  with  the  MR  algorithm  is 
independent  of  image  size. 

The  approach  has  a  number  of  advantages  in  addition  to  the  reduction  in  computational 
cost.  First,  multiresolution  representations  of  the  flow  field  are  available  and,  although  we 
have  not  taken  advantage  of  it  in  this  paper,  the  MR  algorithm  also  allows  for  multireso¬ 
lution  measurements  of  the  optical  flow.  Second,  error  covariance  information  is  available, 
allowing  one  to  assess  the  quality  of  the  estimated  optical  flow.  Finally,  the  MR  algorithm 
is  an  excellent  pre- conditioner  for  algorithms  computing  a  solution  based  on  a  smoothness 
constraint  formulation.  In  addition,  while  we  have  not  pursued  it  here,  the  multiresolution 
philosophy  introduced  here  may  offer  a  promising  approach  to  motion- compensated  image 
sequence  coding.  In  particular,  although  we  used  the  coding  metric  of  reconstruction  error 
as  the  basis  for  the  comparison  of  the  SC  and  MR  approaches  in  Section  3.4,  the  methods 
presented  here  woiild  not  be  the  method  of  choice  for  that  problem.  In  particiileir,  motion- 
compensated  coding  algorithms  designed  specifically  to  minimize  this  criterion  [3,  20,  33] 
will  generally  outperform  the  SC  and  MR  approaches  (which  are  not).  However,  the  com¬ 
putationally  efficient  MR  algorithm  can  be  used  as  an  initial  preconditioning  step  for  such 
coding  algorithms.  In  addition,  one  can  also  imagine  a  second  way  in  which  MR  ideas 
could  be  used  in  this  context.  In  particular,  one  of  the  problems  with  the  SC  and  MR  based 


31 


methods  is  the  differential  form  of  the  brightness  constraint  which,  given  the  discrete  nature 
of  spatial  and  temporfd  sampling,  is  only  valid  for  relatively  small  interframe  displacements; 
In  contrast,  methods  such  as  [3,  20,  33]  use  a  direct  displaced  frame  matching  metric,  which 
is  nothing  but  the  integrated  version  of  the  brightness  constraint.  A  common  approach  to 
dealing  with  larger  displacements  with  the  diflferentM  brightness  constraint  is  to  spatially 
blur  the  image  sequences  -  i.e.  to  consider  lower  resolution  versions  of  the  image  to  estimate 
larger  displacements  [12,  13].  What  this  suggests  is  an  MR  approach  in  which  we  not  only 
have  a  multiresolution  model  for  optical  flow  hut  also  multiresolution  measurements  -  i.e. 
measurements  as  in  (27)  but  for  triples  s  =  (m,  i,  j)  at  several  scales.  The  development  of 
such  an  approach  remains  for  the  future. 

Finally,  in  this  paper  we  have  focused  on  a  particular  image  processing  problem,  the 
computation  of  optical  flow.  However,  we  believe  that  the  multiscale  stochastic  modeling 
approach  can  be  more  generally  useful.  In  particular,  it  may  provide  a  computationally 
attractive  alternative  to  standard  approaches  to  the  broad  class  of  estimation  problems  in 
which  the  underlying  field  to  be  estimated  is  modeled  as  a  Gaussian  Markov  random  field 
or  as  the  solution  of  noise  driven  partial  differential  equations,  or  in  which  a  “smoothness 
constraint”  type  regularization  is  employed.  Viewing  the  multiscale  models  as  an  alternative 
underlying  model  shotdd  lead  to  significant  computational  savings  for  such  problems. 

Acknowledgment;  The  authors  gratefully  acknowledge  the  contributions  of  Dr,  Albert 
Benveniste  who,  among  other  things,  first  suggested  that  the  problem  of  computing  optical 
flow  would  he  well  suited  to  this  approach. 
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A  Non-homogeneous  Tree  Structures 

We  made  the  asstimption  at  the  beginning  of  Section  2.1  that  the  image  lattice  is  square, 
and  that  the  number  of  rows  is  equal  to  a  power  of  two.  The  reason  we  have  done  this 
is  because  of  the  fact  that  the  multiscale  model  described  above  is  defined  on  a  quadtree 
structure.  There  are  at  least  two  ways  to  relax  the  assumption.  First,  we  could  simply  zero 
pad  the  image  lattice  to  make  it  fit  the  quadtree  structure.  In  some  sense,  we  are  changing 
the  image  lattice  to  fit  the  modeling  structure.  A  second,  slightly  more  elegant  approach, 
wotild  be  to  change  the  modeling  structure  to  accommodate  the  lattice.  In  particular,  we 
would  like  to  have  a  structure  which  gives  us  the  proper  number  of  nodes  on  the  finest  level. 
The  quadtree  structure  is  homogeneous  in  the  sense  that  each  parent  has  four  offspring; 
what  we  are  proposing  are  non-homogeneous  tree  structures  in  which  different  parents  may 
have  different  numbers  of  offspring.  For  example,  suppose  one  had  a  6  x  9  lattice.  Figure  27 
illustrates  a  sequence  of  grids  that  one  might  use  to  model  accommodate  this  lattice.  In  the 
first  level,  the  root  node  has  six  offspring,  two  in  the  row  direction  and  three  in  the  column 
direction.  At  the  second  level,  each  node  has  nine  offspring,  three  in  the  row  direction  and 
three  in  the  column  direction.  Thus,  at  the  finest  level  there  is  a  6  x  9  lattice.  This  example 
illustrates  only  one  simple  suggestion.  More  complicated  tree  structures  could  be  derived, 
and  certainly  the  idea  could  be  combined  with  zero  padding. 


B  The  Gauss- Seidel  and  Successive  Over- Relaxation  Algo¬ 
rithms 


One  can  show  using  the  calctilus  of  variations  that  the  solution  to  (4)  satisfies  the  following 
set  of  coupled  partial  differential  equations  [14]: 

=  RE:,^{Et-^VE  -x)  (64) 

V^®2  =  RE^{Et-^VE-x)  (65) 


where: 


E.,  = 
E.,  = 


Et  = 


X  = 
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dzx 
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a2:2 
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E{zx,  Z2,  t) 
E{zi,  Z2,  t) 


dt 

dzi 


E{zi,Z2,t) 


dt 


(66) 

(67) 

(68) 

(69) 


and  where  is  the  Laplacian  operator,  x\  and  X2  are  the  first  €ind  second  components 
of  the  vector  x,  and  R  is  the  parameter  controlling  the  tradeoff  between  the  smoothness 
and  data  dependent  terms  in  (4).  Denote  x{i,j)  =  x{ih,jh,t)  where  h  is  the  grid  spacing. 
Discretizing  (64),  (65)  on  a  uniform  grid  {(i,  j)K  €  {1, ...» Zt},j  €  {1, ...,  Z2}}  leads  to  the 
following  equations  at  each  point: 

—  h  REzyi^EzyXX^i^j  -t"  Ez2^2,i,j  -^t)  (i^O) 
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Figure  27:  Non-homogeneous  tree  structure  for  lattices  which  are  not  square.  The  grid 
structure  is  a  simple  extension  of  the  quadtree  structure  in  that  it  allows  for  varying  numbers 
of  “offspring”  from  each  parent.  The  figure  illustrates  a  hierarchy  of  grids  for  a  6  x  9  lattice. 
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where: 


®2  ~  4a52,i,j  —  h  REz^{Ez^Xi^i^j  +  Ez2^2,i,j  -^t) 


*1  =  +  ®i,tj+i 


®2  =  +  ®2,i+l,i  +  ®2,t,j-l  +  ®2,i,i+l 


(71) 


(72) 

(73) 

(74) 


The  GS  and  SOR  algorithms  separate  the  image  grid  into  two  sets  of  points.  These  are 
generally  referred  to  as  the  Red  points  (*  +  j  is  even)  and  the  Black  points  {i  +  j  is  odd). 
The  Gauss-Seidel  iterations  can  be  derived  by  solving  (70)  and  (71)  for  and  X2,ij: 

GS  Red  Points 


GS  Black  Points 


where: 


—  (®j  —  REz2{Ez2^2,i,3  + 

(75) 

-.n+l 

=  {x^  -  h^REz2iEz,x^+,\  +  Et))/d2,i,3 

(76) 

:"+l 

=  (*?+^  -  h^REz, {Ezzx^tj  +  m/di,i,3 

(77) 

,n+l 

'2,»,j 

=  («?+i  -  h^REz2iEz2^^n  +  Et))/d2,ij 

(78) 

=  4  +  h^REz^ 

(79) 

<3^2, i,i  =  4  +  h^REzz 

(80) 

The  SOR  algorithm  is  very  similar  to  the  GS  algorithm,  except  that  certain  relaxation 
parameters  are  introduced  to  increase  the  convergence  rate.  The  SOR  iterations  are  given 
by: 

SOR  Red  Points 

+«.))/*, (82) 

SOR  Black  Points 

®Ui  =  +  (83) 

®?.tj  =  (l-t«2.i,i)*?,ij  +  a^2,i,i(*?+'-/»^i2^?z.(£?n®"t^  (84) 
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where: 
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(85) 
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Prom  (75),  we  see  that  each  GS  iteration  requires  10  adds  and  4  multiplies  per  pixel 
per  iteration.  Prom  (81),  we  see  that  each  SOR  iteration  requires  12  adds  and  6  multiplies 
per  pixel  per  iteration^®.  Thus,  GS  and  SOR  require  14  eind  18  flops  per  pixel  per  iteration 
respectively. 


C  MR  Algorithm  Complexity  Analysis 

In  this  section  we  analyze  the  computational  complexity  of  the  MR  algorithm.  The  anal¬ 
ysis  applies  to  the  specific  model  given  by  (56)  —  (60).  The  model  is  repeated  here  for 
convenience: 


x{s) 

=  ®(s7)  -j-  (64  )«?(«) 

(89) 

y{^) 

=  C'(s)a!(s)  H- t?(s) 

(90) 

w{s) 

~  Ar{o,i) 

(91) 

v{s) 

~  Ar(0,P(s)) 

(92) 

Xo 

~  //{OypJ) 

(93) 

where  R{s)  =  max(||C'(s)||,  10).  The  analysis  below  takes  into  accotmt  all  floating  point 
adds,  multiplies  and  divides. 

Consider  first  the  update  step  given  by  (42)  -  (45).  P(s|s-|-)  is  initialized  with  pi.  Com¬ 
putation  of  V~^{s)  requires  6  floating  point  operations  (the  inverse  requires  1  divide  since 
F(s)  is  a  scalar  and  the  comparison  required  to  compute  R{s)  is  not  counted).  Computation 
of  K{s)  requires  3  flops.  Computation  of  P(s|s)  requires  7  flops  (Perform  the  C'(a)P(s|s-f-) 
first,  and  use  the  fact  that  P(s|s)  must  be  symmetric).  Initialize  £(s|s-f)  with  zero.  Com¬ 
putation  of  «(s|s)  then  requires  2  flops.  The  update  step  is  required  only  at  the  finest  level, 
since  this  is  the  only  place  we  have  data  for  in  the  optical  flow  problem.  Thus,  the  total 
computation  associated  with  this  step  is  18  X  4^  flops  (I  is  defined  to  be  the  number  of  levels 
in  the  quadtree.  There  are  4^  points  at  the  finest  level.) 

Next,  consider  the  prediction  step,  (46)  -  (48).  Computation  of  Q(saf)  is  negligible 
because  this  parameter  varies  only  as  a  function  of  scale  (level).  Computation  of  P(s|sai) 

All  additions  and  multiplications  which  are  redundant  from  iteration  to  iteration  have  been  ignored.  For 
instance,  1  —  wi.jj  does  not  count  as  an  add  in  (81)  since  one  could  compute  this  once  at  the  start  instead 
of  every  iteration. 
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requires  5  flops  (note  that  F(s)  and  Q(saj)  are  diagonal  multiples  of  the  identity).  Compu¬ 
tation  of  the  predicted  estimate  ®(s|sai)  requires  2  flops.  These  computations  must  be  done 
at  levels  1  through  /.  Thus,  the  total  computation  associated  with  this  step  is  approximately 
7  X  4/3  X  4^  flops. 

Next,  consider  the  merge  step,  (49)  -  (50).  Computation  of  P(s|5-f )  reqtures  44  flops 
(there  are  five  2x2  inverses  requiring  6  flops  apiece,  and  the  computation  of  (1  — 
is  negligible  since  it  only  varies  with  scale.  The  inverses  require  only  6  flops  because  the 
matrices  involved  are  2x2  and  symmetric.)  Computation  of  «(s|s-f)  requires  36  flops.  The 
merge  step  must  be  done  at  levels  0  through  /  —  1.  Thus,  the  total  computation  associated 
with  this  step  is  80  x  1/3  X  4^  flops. 

Finally,  consider  the  steps  in  the  downward  sweep,  (53)  —  (55).  Computation  of  J{s) 
requires  12  flops  (the  matrix  P(s7|s)  has  already  been  inverted  in  (50),  F{s)  is  a  multiple 
of  the  identity  and  J(s)  is  symmetric.)  Computation  of  P»(s)  is  not  required,  unless  one 
is  explicitly  interested  in  the  error  covariance  of  the  smoothed  estimate.  Computation  of 
*,(5)  requires  10  flops.  The  smoothing  step  must  be  done  at  levels  1  through  /.  Thus,  the 
total  computation  associated  with  this  step  is  22  x  4^  flops. 

We  can  now  add  up  aU  of  the  computations  associated  with  the  MR  algorithm.  There 
are  4*  pixels  in  the  problem  domain,  and  thus  the  algorithm  requires  18  4-28/3-1-  80/3  -f 
22  =  76  flops  per  pixel.  We  note  that  this  is  a  lower  bound  on  the  number  of  flops  per 
pixel  in  any  implementation  of  the  algorithm  and  that  the  implementation  with  the  lowest 
number  of  flops  per  pixel  may  not  be  the  best.  The  reason  is  simply  that  there  may  not  be 
enough  memory  available  to  keep  aU  intermediate  calculations  around  (such  as  the  inverses 
computed  in  (50)  and  reused  in  (55)).  We  compute  the  complexity  of  the  GS  and  SOR 
algorithms  in  the  same  way  (i.e.  aU  intermediate  results  are  assumed  to  be  available),  and 
thus  the  computational  comparison  we  make  between  these  algorithms  is  based  on  optimal 
(in  terms  of  the  number  of  flops)  implementations.  Suboptimal  implementation  of  the  MR 
algorithm  will  lower  its  computational  advantage,  but  any  reasonable  implementation  (for 
instance  one  which  saves  just  ®(s|5),  P(a|5)  and  the  measurement  data)  will  still  provide  a 
significant  savings  over  the  SOR  and  GS  algorithms. 


37 


References 

[1]  J.  Aggarwal  ano  N.  Nandhakumar,  “On  the  computation  of  motion  from  se¬ 
quences  of  images  -  a  review.”  Proc.  IE£E,  76:917-935, 1988. 

[2]  P.  Anandan,  “A  computationad  framework  and  an  algorithm  for  the  measurement  of 
visual  motion,”  International  Journal  of  Computer  Vision,  2:283-310, 1989. 

[3]  N.  Baaziz  and  C.  Labit,  “Multigrid  motion  estimation  on  pyramidal  representations 
for  image  sequence  coding,”  IBISA  internal  publication  No.  572,  February  1991. 

[4]  M.  Bbrtero,  T.  Poggio  and  V.  Torre,  “Hi-posed  problems  in  early  vision,”  Proc. 
IEEE,  76:869-889, 1988. 

[5]  P.  Burt  and  E.  Adelson,  “The  Laplacian  Pyramid  as  a  compact  image  code,”  IEEE 
Trans.  Comm.y  31:482-540, 1983. 

[6]  T.  Chin,  “Dynamic  Estimation  in  Computational  Vision,”  MIT  Dept,  of  EECS  Ph.D. 
Thesis,  Oct.  1991. 

[7]  K.  C.  Chou,  A.  S.  Willsky,  A.  Benveniste  and  M.  Basseville,  “Recursive  and 
Iterative  estimation  algorithms  for  multiresolution  stochastic  processes,”  Proc.  of  the 
IEEE  CDC,  Dec.  1989. 

[8]  K.  C.  Chou,  A  stochastic  modeling  approach  to  multiscale  signal  processing,  MIT 
Dept,  of  EECS  Ph.D.  Thesis,  May  1991. 

[9]  K.  C.  Chou,  A.  S.  Willsky  and  A.  Benveniste,  “Multiscale  Recursive  Estimation, 
Data  Fusion  and  Regularization,”  MIT  LIDS  Report  ^  2085,  submitted  for  publication 
to  IEEE  Trans,  on  Automatic  Control. 

[10]  K.  C.  Chou,  A.  S.  Willsky  and  R*  Nikoukhah,  “Multiscale  Systems,  Kahnan 
Filters  and  Riccati  Equations,”  submitted  for  publication  to  IEEE  Trans,  on  Automatic 
Control. 

[11]  S.  C.  Clippingdale  and  R.  G.  Wilson,  “Least  Squares  Image  Estimation  on  a  Mul¬ 
tiresolution  Pyramid,”  Proceedings  of  the  1989  International  Conference  on  Acoustics, 
Speech  and  Signal  Processing. 

[12]  W.  Enkelmann,  “Investigations  of  multigrid  algorithms  for  the  estimation  of  optical 
flow  fields  in  image  sequences,”  Computer  Vision,  Graphics  and  Image  Proc.,  43:150- 
177,  1988. 

[13]  Heitz,  F.  and  Bouthemy,  P.,  Multimodal  estimation  of  discontinuous  optical  flow 
using  Markov  random  fields,  INRIA  Report  1367,  Jem.  1991,  eilso  submitted  to  IEEE 
Trans,  on  Pattern  Analysis  and  Machine  Intelligence. 

[14]  B.  K.  P.  Horn  and  B.  Schunck,  “Determining  Optical  Flow,”  Artificial  Intelligence, 
17:185-203, 1981. 


38 


[15]  J.  Kearney,  W.  Thompson  and  D.  Boley,  “Optical  Flow  Estimation:  An  Error 
Analysis  of  Gradient  Based  Methods  with  Local  Optimization,”  IEEE  Trans.  PAMI, 
2:229-244, 1987. 

[16]  G.-C  J.  Kuo,  B.  C.  Levy  and  B.  R.  Musicus,  “A  local  relaxation  method  for 
solving  elliptic  PDE’s  on  mesh  coimected  arrays,”  SIAM  J.  Sci.  Stat.  Comput.,  8:550- 
573,  1987. 

[17]  S.  Mallat,  “Multi-frequency  Channel  Decomposition  of  Images  and  Wavelet  Models,” 
IEEE  Trans.  ASSP,  37:2091-2110, 1989. 

[18]  B.  Mandelbrot  and  H.  Van  Ness,  “Fractional  Brownian  Motions,  Fractional  Noises 
and  Applications,”  SIAM  Review,  10:422-436,  Oct.  1968. 

[19]  D.  Murray  and  B.  Buxton,  “Scene  Segmentation  From  Visual  Motion  Using  Global 
Optimization,”  IEEE  PAMI  2:220-228, 1987. 

[20]  A.  Netravali  and  J.  Robbins,  “Motion- Compensated  Television  Coding:  Part  1,” 
Bell  System  Technical  Journal  58:631  -  670,  1979. 

[21]  J.  Prince  and  E.  McVeigh,  “Motion  estimation  from  tagged  MR  image  sequences.” 
Johns  Hopkins  University  Dept,  of  Elec,  and  Computer  Engineering  Technical  Report 
JHU/ECE  91-05,  also  submitted  to  IEEE  Transactions  on  Medical  Imaging,  1991. 

[22]  D.  Raviv,  “A  Quantitative  Approach  to  Camera  Fixation,”  in  Proc.  IEEE  Conf.  on 
Computer  Vision  and  Pattern  Recognition,  Maui,  Hawaii,  1991. 

[23]  A.  Rougee,  B.  Levy  and  A.  S.  Willsky,  “An  estimation-based  approach  to  the 
reconstruction  of  optical  flow,”  Laboratory  for  Information  and  Decision  Systems  Tech¬ 
nical  Report,  No.  1663,  MIT,  April  1987. 

[24]  A.  Rouged,  B.  Levy  and  A.  S.  Willsky,  “Reconstruction  of  two  dimensional 
velocity  fields  as  a  linear  estimation  problem,”  in  Proceedings  1st  Int.  Conf.  Computer 
Vision  (London,  England)  pp.  646-650,  June  1987. 

[25]  A.  Shio  and  j.  Sklansky,  “Segmentation  of  People  in  Motion,”  in  Proc.  IEEE 
Workshop  on  Visual  Motion,  Princeton  NJ,  Oct.  7-9,  1991,  pp.  325  -  332. 

[26]  E.  SiMONCELLi  AND  E.  Adelson  AND  D.  Hbeger,  “Probability  Distributions  of 
Optical  Flow,”  IEEE  Conference  on  Computer  Vision  and  Pattern  Recognition,  Maui 
Hawmi,  June  1991. 

[27]  G.  Strang,  Introduction  to  A f plied  Mathematics.  WeUesley-Ceimbridge  Press,  Welles¬ 
ley  MA,  1986. 

[28]  G.  Strang,  Linear  Algebra  and  its  Applications.  Harcourt-Brace-Jovanovich,  1988. 

[29]  V ;  Sund  ARBS WAREN ,  “Egomotion from  Global  Flow  Field  Data,”  in  Proc.  IEEE  Work¬ 
shop  on  Visual  Motion,  Princeton  NJ,  Oct.  7-9,  1991,  pp.  140-5. 


39 


[30]  R.  SZELISKI,  Bayesian  Modeling  of  Uncertainty  in  Low-level  Vision.  Kluwer  Academic, 
1989. 

[31]  D.  Terzopoulos,  “Image  analysis  using  multigrid  relaxation  methods,”  IEEE  Trans, 
on  PAMI,  8:129-139,  1986. 

[32]  H.  L.  VanTrbes,  Detection,  Estimation  and  Modulation  Theory:  Part  1.  Wiley,  New 
York,  1968. 

[33]  D.  Walker  and  K.  Rao,  “Improved  pel-recursive  motion  compensation,”  IEEE 
Trans,  on  Communications,  32:1128-1134, 1984. 

[34]  G.  WoRNELL,  “A  Karhtmen-Loeve  like  expansion  for  1/f  processes,”  IEEE  Trams,  on 
Inf.  Theory,  36:859-861,  1990. 


40 


